function [out coid CC] = find_in_ct(path, panelHandle)
% Function search for Electrodes points in CT image
% path - path to CT analyze image header file (.hdr)
% out.path - same path
% out.EL - matrix of provisional electrodes point coordinates (voxel space) in CT

% nfo = analyze75info(path);
% X1 = analyze75read(nfo);
% Y1 = mat2gray(X1);
% P = im2uint8(Y1);
% clear Y1 X1
stack = get(panelHandle,'userdata');
P = stack(7).vol;
%P = mat2gray(im2uint8(P));


s = size(P);
%% Find intracerebral space
H = logical(zeros(s));

se = strel('disk', 10);
se1 = strel('disk', 5);
seD = strel('disk',30);
P = rescale_ct(P);
bone_value = 170;
electrodes_value = 230;

for i=1:s(3)
    a = squeeze(P(:,:,i));
    r = a;
    a(a>=bone_value) = 255;
    a(a~=255) = 0;
    b = im2bw(a);
    
    
    c = imfill(b, 'holes');
    c = imclose(c, se1);
    c = imfill(c, 'holes');
    statsi = regionprops(c, 'Area');
    A = [statsi.Area];
  if(isempty(A))
        H(:,:,i) = c;
        continue;
    end
    A = sort(A);

    c = bwareaopen(c, floor(A(end)*0.3));
    c = imdilate(c, se);

    
    c = imfill(c, 'holes');
    c = imerode(c, se);
    
    stats = regionprops(c, 'Area','Perimeter');
    A = [stats.Area];
    [A ix] = sort(A);
    Per = [stats.Perimeter];
    Per = Per(ix);
    u = 4*pi*A(end)/Per(end)^2;
    tempA = A(end);
    if u < 0.5
        
        
        BB = regionprops(c, 'BoundingBox');
        CH = regionprops(c, 'ConvexImage');
        bb = BB.BoundingBox;
        ch = CH.ConvexImage;
        
        r = c;
        r(floor(bb(2))+1:floor(bb(2))+bb(4), floor(bb(1))+1:floor(bb(1))+bb(3)) = CH(1).ConvexImage;
        
        r1 = r - c;
        L = bwlabel(r1);
        stats = regionprops(L, 'Area');
        A = [stats.Area];
        [A] = sort(A);
        
        if(~isempty(A))
            if A(end) > tempA
                r2 = bwareaopen(r1, A(end)-1);
                
                c = r2+c;
            end
        end
    end
    
    
    b = imopen(c, seD);
    H(:,:,i) = b;
end
% 
% coid = regionprops(H, 'Centroid');
% coid = coid.Centroid;

coid = regionprops(H, 'Centroid', 'Area');
[v, i] = max([coid.Area]);
coid = coid(i).Centroid;
%% Do segmentation
P1 = P.*uint8(H);

s = size(P1);

H1 = logical(zeros(s));

for i=1:s(3)
    
    a = squeeze(P1(:,:,i));

    a(a>=electrodes_value) = 255;
    a(a~=255) = 0;
    b = im2bw(a);

    H1(:,:,i) = b;

end

H1 = imdilate(H1, ones(5,5,5));

%% Do Correlation

P2 = P.*uint8(H1);

I2 = double(P2);
clear P2
s = size(I2);
I1 = zeros(s);
load strelem

% pp = 4;
I1(round(end/2)-1:round(end/2)+3,round(end/2)-1:round(end/2)+3,round(end/2)-1:round(end/2)+3) = G;
pp = 0;
n = 3;
% I1(pp+round(end/2)-n:pp+round(end/2)+n,pp+round(end/2)-n:pp+round(end/2)+n,pp+round(end/2)-n:pp+round(end/2)+n) = G7;

CC = ifftn(fftn(I2).*conj(fftn(I1)));
clear I1 I2
CC = fftshift(CC);
CC = mat2gray(CC);

fgm = imregionalmax(CC.*double(H1));

%% Find electrode points
s = size(fgm);

X = 0;
Y = 0;
Z = 0;

for i=1:s(3)

   [x y] = find(fgm(:,:,i) == 1);

   if ~isempty(x)
       z = ones(size(x))*i;
       X = [X; x];
       Y = [Y; y];
       Z = [Z; z];
   end
end


X = X(2:end);
Y = Y(2:end);
Z = Z(2:end);

%% Save result
out.path = path;
out.El = [X(:) Y(:) Z(:)];